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Abstract.  We  consider  a  class  of  probability  measure  dependent  dynamical  systems  which  arise 
in  the  study  of  multiscale  phenomena  in  diverse  fields  such  as  immunological  population  dynam¬ 
ics,  viscoelasticity  of  polymers  and  rubber,  and  polarization  in  dielectric  materials.  We  develop  an 
inverse  problem  framework  for  studying  systems  with  distributed  temporal  delays.  In  particular, 
we  establish  conditions  for  existence  and  uniqueness  of  the  forward  problem  and  well-posedness 
(including  method  stability  under  numerical  approximations)  for  the  inverse  problem  of  estimating 
the  probability  measures.  We  show  that  a  motivating  class  of  models  of  HIV  infection  dynamics 
satisfies  all  the  conditions  of  our  framework,  thereby  providing  a  theoretical  foundation  for  inverse 
problem  computations  with  these  models. 


1.  Introduction 

In  this  paper,  we  present  a  general  theoretical  framework  including  implementable  approxima¬ 
tion  ideas  for  inverse  problems  involving  measure  dependent  dynamical  systems.  In  particular,  we 
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treat  inverse  problems  for  systems  wherein  the  dynamics  are  given  for  the  expected  values  of  some 
of  the  states 

x(t,P)  =  gP[x(t,q)\  =  /  x(t,q)dP(q) , 

Jq 

with  respect  to  a  probability  measure  P  defined  on  the  Borel  measurable  subsets  of  an  admissible 
parameter  set  Q.  Here,  for  example,  the  expected  state  dynamics  might  be  given  by  a  differential 
equation 


*0)  =f(t,x(t),P), 

where  the  right  side  is  dependent  on  the  measure  P.  In  general,  the  function  /  can  also  represent 
a  delay  or  partial  differential  equation.  Whatever  the  form  of  the  system,  the  individual  dynamics 
(for  each  q  e  Q)  for  x(t.q)  are  not  available  in  the  cases  of  interest  to  us  here. 

We  specifically  focus  on  a  framework  for  systems  of  measure  dependent  delay  differential  equa¬ 
tions  that  are  motivated  by  examples  arising  in  cellular-level  modeling  of  HIV  pathogenesis.  In 
this  case,  the  probability  distributions  characterize  delays  inherent  in  cellular  pathways  for  virus 
production.  These  delays  represent  lags  between  acute  infection  of  cells  and  the  initiation  of  vi¬ 
ral  production  as  well  as  between  productive  acute  infection  and  chronic  infection  (e.g.,  see  the 
discussions  in  [19,  22,  25,  29,  30,  31,  32]).  The  particular  model  that  motivates  our  mathematical 
efforts  here  was  derived  in  [6]  to  describe  infection  dynamics  in  an  in  vitro  experiment  and  is  given 
by  the  system  of  integro-differential  equations 

r° 

(1.1)  V(t)  =  -cV(t)+nA  A(t  +  x)dPi{x)+ncC{t)-pV{t)T{t ) 

r° 

(1.2)  A(t)  =  (rv-8A-8X(t))A(t)-yl  A(t  +  x)dP2{x)  +  pV(t)T(t) 

fO 

(1.3)  C(t)  =  (rv-8c-8X(t))C(t)  +y J  A(t  +  x)dP2{x) 

(1.4)  T(t)  -  (ru  -  K  -  6 X(t)  - pV(t ))  T(t)  +S. 
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The  compartment  variables  consist  of  the  virus  V  and  the  acutely  infected,  chronically  infected, 
and  uninfected  T-cells,  A,  C,  and  T,  respectively,  while  X  represents  the  total  cell  population 
(A  +  C  +  T).  The  probability  measures  P\  and  Pi  correspond,  respectively,  to  the  distributions  of 
the  two  pathogenic  delays  discussed  above.  In  the  derivation  of  this  model,  we  ascribed  no  subclass 
variability  to  A  and  T  and  consequently  only  the  variables  C  and  V  are  actually  expected  values 
(with  respect  to  Pi  and  Pi,  respectively).  All  constitutive  parameters  (c,  ua,  nc,  p,  rv,  8a,  8,  y,  8c, 
ru,  8U )  and  details  of  the  model  as  well  as  its  derivation  are  fully  described  in  [6]. 

The  example  motivating  our  efforts  here  is  merely  one  of  several  important  examples  of  measure 
dependent  dynamics  that  have  arisen  in  recent  applications.  Realistic  models  in  viscoelasticity  of 
polymers  and  rubber  as  well  as  polarization  in  dielectric  materials  share  certain  features  with  our 
formulation.  Specifically,  recent  studies  [12,  13]  of  molecular-based  stick-slip  reptation  models 
for  heterogeneous  viscoelastic  polymer  chain  materials  involve  systems  of  the  form 


(1.5) 


d2u(t,x) 
Y  dt2 


d_ 

dx 


c(t,x;P)  —  F(t) , 


where  u  is  the  tensile  displacement,  o  is  the  measure  dependent  stress 


(1.6) 


o(t,x;P)  =£e(e(f),e(f))+v 


Ei(t,x\z)dP(%) , 


with  strain  £  =  ^  and  “internal”  strain  £i  defined  by 


(1.7) 


0£l 

~di 


(t,x;x)  +  -£i(t,x;x) 

T 


e(t,x)h(£(t,x)) . 


In  [14]  it  is  shown  that  similar  models  (with  different  nonlinearities  in  (1.6),  (1.7))  are  important 
if  one  replaces  Fung  kernels  [21]  with  equations  for  distributed  molecular  mechanisms  in  describ¬ 
ing  shear  response  in  biotissue.  In  another  important  application  [8],  the  systems  are  the  usual 
Maxwell’s  equations  for  the  electromagnetic  fields  E  and  H  in  a  heterogeneous  dielectric  and  are 
given  by 

Vx£  =  VD  =  0 

VxH  =  +j  S7H  =  0 
D  =  £,  E  +  P , 
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with  the  exception  that  the  probability  measure  (P)  dependent  macroscopic  polarization  P  is  given 
by 

P(t)  =  [  [pi(p,x)+p2{t;d)\dP{x,d). 

J'Tx'D 

In  this  case,  the  microscopic  orientational  (Debye)  polarization  pi(t',x)  is  defined  by 

1 

P 1  +  ~Pi  =  e£, 

T 

and  the  microscopic  electronic  (Lorentz)  polarization  p2(f,d)  is  defined  by 

mp2  +  cp2  +  kp2  =  eE , 

with  d  =  — .  In  the  equations  for  both  p\  and  p2,  the  parameters  x  and  d  represent  relaxation 
parameters  which  may  vary  over  the  admissible  set  CT  x  P>  according  to  some  unknown  but  sought 
after  probability  measure  P  —  P(x.d). 

In  each  of  these  examples,  one  seeks  to  characterize  the  material  behavior  to  perturbations  by 
finding  a  measure  P*  that  provides  the  “best”  mathematical  system  response  when  compared  to 
observations  of  the  physical  system. 

Finally,  a  more  classical  use  of  measure  dependent  dynamical  systems  can  be  found  in  the 
literature  on  “relaxed”  or  “sliding  regime”  control  systems  (e.g.,  [28,  35,  36]  as  well  as  in  Preisach 
models  for  hysteresis  in  materials  (e.g.,  see  [10,  11,  24,  27,  34]  and  the  references  therein). 

In  this  paper,  we  formulate  an  inverse  problem  framework  for  estimating  the  measures  (us¬ 
ing  system  observations)  in  systems  of  probability  measure -dependent  nonlinear  delay  differential 
equations.  To  do  this,  we  employ  the  Prohorov  metric  (equivalent  to  weak  convergence  of  mea¬ 
sures)  in  a  functional  analytic  setting.  We  also  show  how  to  develop  an  approximation  approach 
that  is  readily  implemented  to  obtain  computational  methods.  Finally,  we  illustrate  that  the  HIV 
systems  in  [6]  are  concrete  realizations  of  systems  satisfying  the  hypotheses  in  our  framework.  We 
thus  provide  a  theoretical  foundation  for  the  computational  methods  used  with  experimental  data 
in  [6]  as  well  as  provide  a  general  framework  for  a  class  of  systems  arising  in  other  applications. 
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2.  Theoretical  Framework  for  General  Delay  System 
2.1.  Theoretical  Framework.  We  begin  by  considering  an  n-dimcnsional  generalized  delay  sys- 


tern: 

(2.1) 

x(t) 

=  I,.  f’l.P; - P„,,l  .1^0 

(2.2) 

(x(0),x0) 

=  (<J>(0),K), 

for  the  piecewise  absolutely  continuous  (PAC)  initial  condition  <t>,  defined  on  the  hysteretic  do¬ 
main  Q  —  [— r,0]  where  r  >  0.  Following  the  standard  practice  in  the  delay  differential  equation 
literature,  we  let  the  notation  xt  represent  the  function  0  i— >•  x(t  +  0),  —  r  <  0  <  0.  Furthermore, 
we  let  §  be  the  class  of  all  Borel  subsets  of  Q  (the  Borel  o-algebra),  and  £?(Q)  be  the  space  of 
probability  measures  on  (<2,§).  In  order  to  make  this  framework  applicable  to  realistic  problems, 
we  restrict  the  space  of  admissible  probability  distributions  ad{Q )  to  PAC  functions  with  a  finite 
number  of  saltations,  where  P,  e  & ad(Q ),  i  —  1.2.....  np.  To  establish  the  theoretical  framework 
necessary  to  estimate  a  finite  number  np  of  probability  distributions  P,,  i  =  1 , . . . ,  np,  in  (2.1),  we 
define  the  set 

nad  =  f[  &ad(Q). 

i=  1 

Thus  n  E  nad  means  that  n  =  (Pi , . . . ,  PHp ),  for  P] , . . . ,  Pl2p  e  & ad  ( Q )  • 

In  the  following  discussion,  we  make  use  of  a  construct  from  advanced  probability  theory,  the 
Prohorov  metric  (denoted  by  p).  As  we  note  below,  convergence  in  the  Prohorov  metric  is  equiva¬ 
lent  to  “weak  convergence  of  measures”  (which  is  actually  weak*  convergence  when  considering 
&  C  ho* ,  where  is  the  space  of  continuous  functions  on  Q  with  the  max  norm).  We  recall  the 
definition  of  the  Prohorov  metric  on  tP(Q)  for  any  arbitrary  complete  metric  space  Q  with  metric 
dist.  For  any  closed  F  C  Q  and  £  >  0,  we  define  the  £-neighborhood  of  F  by 

F£  =  {q  G  Q  :  dist(q,q)  <£,  for  some  q  £  Fj  . 
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The  Prohorov  metric  is  then  the  mapping  p  :  &(Q)  x  dP(Q)  — >  M+  defined  by 

p {Pi, Pi)  =  inf{£  >  0  :  P\(F)  <  P2(FE)  +  £,  F  closed,  F  C  <2}  . 

It  is  known  (see  [16,  17,  20,  26,  33])  that  convergence  Pn  — >  P  in  the  Prohorov  metric  is  equivalent 
to  the  statement 

[  §dPn  — >  [  §dP  for  all  (f)  £  %?(Q) , 

Jq  Jq 

which  is  convergence  in  expectation  or  distribution,  also  called  weak  convergence  of  measures. 
We  define  a  topology  for  Ilad  by  extending  the  Prohorov  metric  such  that  for  71,  ft  6  11^,  we  let 

np 

Pn(rc,Tt)  =  £p  {Pu Pi), 

i=  1 

where  n  =  (Pi,. . . ,  P,lp)  and  Tt~(P\,.. . .  P,lp )  •  Given  our  data  dj  e  Wld,  nd  <  n,  sampled  at  discrete 
times  tj,  i  =  1,2,...  ,nt  and  observation  matrix  C  £  R'vx'\  the  goal  of  the  inverse  problem  is  to  find 
a  solution  to 

nt 

(2.3)  min  /(7t;d)  =  min  V  \Cx(tj.K)  —dj\2  , 

where  y  is  the  solution  to  (2.1)-(2.2)  corresponding  to  n  and  d  =  {d\,  di_,  ■ . .  d,y}.  Note  that  in 
general,  J  need  not  have  a  unique  minimizer,  and  thus  the  corresponding  solution  (denoted  by 
Il*(d))  could  be  a  set  of  probability  distribution  functions.  In  this  case,  we  would  then  define  the 
distance  between  two  of  these  sets  II*  (d)  and  n*(d)  (for  data  d,d)  to  be 

dH( n*(d),n*(S))  =  inf  {pn(7C,7t)  :  n  £  Il*(d),fc  £  n*(d)}  , 

which  is  the  well-known  Hausdorff  distance  (see  [23]). 

To  establish  the  well-posedness  of  (2.3)  within  our  framework,  we  now  examine  the  forward 


and  inverse  problems. 
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2.2.  Forward  Problem.  To  prove  the  well-posedness  of  the  forward  problem,  we  need  to  estab¬ 
lish  conditions  for  the  existence  and  uniqueness  of  a  solution  to  (2.1)-(2.2)  as  well  as  its  continuous 
dependence  upon  n  in  the  pn  metric.  We  first  give  conditions  on  the  function  g  in  (2.1). 

Condition  2.1.  Let  (rj, (|), tc) ,  (fj,  (j). 7t)  e  M'!  x  PAC(0;R'!)  x  Had-  We  require  that  the  function 
g  :  1RW  x  PAC(2;K”)  x  Uad  — »  M",  satisfy 

where  K[  >  0,  j  •  |,  •  |  arc  the  Euclidean  norm  in  R"  and  the  max  norm  on  PAC((2;RR),  respectively, 
and  the  last  term  on  the  right  side  <^((]);7t,7t)  is  some  function  such  that  for  each  (f)  €  PAC((2;R'7), 
^  0  as  pn(7t,fc)  — >  0  and  <  K\§ |oo. 


Following  the  standard  Picard  iteration  arguments  for  proving  the  existence  and  uniqueness  of 
the  solution  to  an  ordinary  differential  equation  on  a  finite  interval  I  —  [0,  tj] ,  we  begin  by  rewriting 
(2.1)-(2.2)  as 

(  ^>(0)  +  /og(ac(s),xs,7t)ds  for  tel 
(2.4)  x{t)  =  l 

{  ®(t)  for  teQ, 

with  initial  condition  <f>.  In  the  subsequent  proof,  we  will  make  use  of  the  following  definition. 


Definition  2.2.  Let  successive  approximations  to  the  solution  of  (2.4)  on  [—r,tf]  be  defined  for 
j  =  0,1,2,...,  as 


(2.5) 


yj+\t) 


d>(0)  for  tel 
d>(t)  for  teQ 

<p(0)  +  fog(yJ(s),yi,7t)ds  for  tel 
d>(?)  for  teQ, 


ad- 


with  initial  condition  <t>  and  7t  e  fl, 
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Theorem  2.3.  Consider  the  system  of  equations  (2.1)  with  initial  condition  0  E  PAC(Q:W')  and  t 
in  a  finite  interx’al  I  =  [0  ,tf],  0  <tf  <°°.  If  the  function  g  satisfies  Condition  2.1,  then  there  exists 
a  unique  solution  to  (2.1)  on  I. 

Proof.  The  general  idea  of  our  argument,  which  is  quite  standard,  is  to  show  that  the  successive 
approximations  defined  in  (2.5)  converge  to  a  unique  solution  of  (2.1). 

Let  the  residual  function  of  two  functions  xt,  xt  be  defined  as 

(2.6)  e(t;xt,xt )  =  \x(t)  ^ x(t)\  +  \xt  - x^  , 
for  xt,xt  E  PAC  (<2;MW),  t  El. 

If  we  consider  the  residual  for  the  functions  yJ+l  and  yf  as  defined  in  (2.5)  we  find  that  for  t  El 
and  j  >  0 

e{t;yj+l,y{)  =  {g(/ (5), yi,n) -g(yJ~l(s):yJ~l,Ti)}ds 

+  [  {g(yJ(s),yt,Tl)-g(yJ-\s),yt-\Tl)}ds 
J  0  00 

<  KL  [  e(s;yJs,yJ~l)ds+  KL  [  e(s;yj,yJ~l)ds 

Jo  Jo 

and  thus 

(2.7)  e(t;y{+l ,y{)  <  2KL  [  e(s;yJs,yJ~x)ds . 

Jo 

Note  that  the  case  for  j  —  0  (with  t  El)  is  special: 

=  |/(0-y0(0|  + |yf-y?L 

=  f  g(y\S),y0s,K)ds  +  [t+  g(y0(s),y°s,n)ds 
Jo  Jo 

<  2  J  {^L|e(5,y°,0)|}  ds 

<  2KL  f  (\<P(s)\  +  \<PUds 

Jo 

<  2ZL(|0(O)|  +  |0|Jt, 


INVERSE  PROBLEMS  FOR  A  CLASS  OF  MEASURE  DEPENDENT  DYNAMICAL  SYSTEMS 


9 


and  thus 


e(t-yj.yi  )  <  KGt , 


where 


Kg  =  2^  (1^(0)  |  +  |d>| 


We  claim  that  from  (2.7)  and  the  j  =  0  case,  we  have 


e(t-,y}+1  ,y{)  < 


Kg  (2 KLt)j+ 
2 Kl  (j  +  1)! 


■,  tel. 


Clearly,  this  is  true  for  j  —  0,  and  the  general  case  follows  easily  from  induction  using  (2.7).  Using 
the  estimate  (2.8),  we  can  then  infer  that 

f  j+l  a  ^  (2Ku)V+V 

,y,)  - 

<  E^z2KLt. 

-  2  Kl 

Thus,  by  the  comparison  test,  £J=0e(t;y/+1,y/)  converges  uniformly  for  tel,  which  proves  that 
{yj(t)}  converges  uniformly  for  all  tel.  Denote  lim/_^ooy-7(t)  asy(f).  Since  then’s  are  continuous 
and  converge  uniformly  to  y,  we  see  that  y  is  both  continuous  on  /  and  satisfies  (2.1)  by  taking  limits 
in  (2.5).  Note  that  this  also  yields  y  absolutely  continuous  on  I. 

To  prove  the  uniqueness  of  our  solution,  suppose  we  have  two  distinct  solutions  {y,y}  e  ^(7) 
to  (2.1).  Using  the  same  arguments  as  in  establishing  (2.7),  we  have 


e(t;ytJt)  <  /  {g(y(s)++)-g(X5)++))}^ 

Jo 


+  /  {g(y(s)+,TC)-g(:v(s)++))}^ 

Jo 


<  2 Kl  [  e(s;ys,ys)ds . 
Jo 


Thus  by  Gronwall’s  inequality  we  have  that 


|y(t)-y(OI  +  |vr-yf|oo  <0  for  tel. 
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Thus  y(t)  =  y(t)  for  16  /  and  also  for  t  £  Q  since  both  solutions  satisfy  the  same  initial  condition. 

We  have  therefore  now  proven  that  there  exists  a  unique  solution  (for  t  £  [— r,//])  10  (2.4)  and 
thus  to  (2.1),  and  that  the  solution  is  in  fact  absolutely  continuous  on  I.  □ 

Remark  2.4.  Existence  and  uniqueness  can  be  established  on  any  finite  interval  I  under  somewhat 
weaker  conditions  than  the  global  Lipschitz  requirement  in  Condition  2.1.  In  particular,  under 
conditions  of  local  Lipschitz  plus  affine  growth  domination  at  infinity,  one  can  also  establish  the 
existence  and  uniqueness  results  of  Theorem  2.3.  For  examples  and  details  of  such  arguments  in 
the  case  of  general  nonlinear  systems,  see  [1,  2,  15]. 

In  order  to  consider  calculating  a  solution  to  the  inverse  problem  (2.3),  we  must  establish  the 
well-posedness  of  this  optimization  problem,  which  involves  both  proving  that  the  solution  to  (2.1) 
is  continuous  in  n  and  that  there  exists  a  solution  to  (2.3).  Following  the  development  in  [7],  with 
inspiration  from  [4,  5],  we  will  do  so  by  examining  both  problem  stability  and  method  stability. 

We  say  that  the  forward  problem  is  well-posed  if  the  unique  solution  to  the  model  (in  our  case 
x )  is  continuously  dependent  upon  the  measures  K.  Let  us  fix  t  £  I  and  consider  the  continuity 
of  a  solution  to  (2.1)-(2.2)  with  respect  to  n.  Thus  we  interpret  the  solution  x  as  the  mapping 
x(t,  •)  :  Ilad  — ■>  R",  parameterized  by  the  time  t. 

Lemma  2.5.  If  l  £  /  =  [0,  tf\,  Q  =  [— r,  0],  tl  £  Ylad>  and  the  right  side  g  of  (2.1)  satisfies  Condition 
2.1,  then  the  unique  solution  to  (2.1)-(2.2)  is  point-wise  continuous  at  K  £  Ylad. 

Proof.  Since  p  is  a  metric  topology,  it  suffices  to  argue  that  x(t,ltj)  — >  x(/,7t)  for  any  sequence 
{Hi}f=l  where  as  i  —>  °°,  tij  — >  n  in  pn-  From  the  definition  of  solutions,  we  find 

|*(f,7C,-)-*(f,7C)|  <  /  )g(x(s,7l;),x,(7t*),7q-) -g(*(s,7l),*lS(7C),7C)|  ds . 

Jo 

Since  g  satisfies  Condition  2.1,  we  have  that 

(2.9)  \x(t,Ki)-x(t,K)\<  f  KLe(s;xs(Tli),xs(Tl))ds  +  &(Tli:Tl), 

Jo 
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where  e  is  the  residual  as  defined  in  (2.6)  and  it)  =  fff  Next,  consider 

rt+- 


|^(7Ci)-^(7C)|00  < 


1 0 


| g(x(s,Tti),xs(%i),Tti)  -g(x(s,n)1xs(n),Ti)\ds 


<  /  \8{x(s,ni)1xs(TZi)1ni)-g(x(s,n)1xs(n),n)\ds, 

Jo 

which  has  the  same  bound  as  in  (2.9).  By  combining  these  bounds  we  find  that 


e(t;xt(ni),xt(n))  <  2 KL  J I  e(s-,xs(ni),xs(%))ds 


+2«^(7Cf,7C), 


and  an  application  of  Gronwall’s  inequality  yields 

e(t;xt(ni),xt(n))  <  2^(ni,n)efo2KLds 
<  2‘7(ni,n)e2KLtf . 

From  this  it  immediately  follows  that  \x(t,Ki)  —  x(t.7l)\  — >•  0  as  i  — >  °o  and  7t,  — ^  tt  in  the  Prohorov 
metric.  Therefore,  we  have  pointwise  continuity  of  the  solution  x  (and  thus  J)  with  respect  to  the 
optimization  variable  n  of  interest  in  (2.3).  □ 


2.3.  Inverse  Problem  Stability  and  Method  Stability.  As  mentioned  at  the  beginning  of  the  last 
section,  in  order  to  fully  justify  our  claim  regarding  the  well  posedness  of  the  inverse  problem,  we 
need  to  examine  questions  concerning  the  existence  of  a  solution  to  (2.3)  as  well  as  the  dependence 
of  those  solutions  upon  given  data. 


Theorem  2.6.  There  exists  a  solution  to  the  inverse  problem  as  described  in  (2.3). 

Proof.  From  results  in  [17],  we  know  that  if  Q  is  compact,  (^ad(Q)-P)  is  compact  and  thus  by 
extension,  (II ad,  pn)  is  also  compact.  It  is  well  known  that  a  continuous  function  on  a  compact  set 
in  a  metric  space  attains  a  maximum  and  a  minimum.  By  Lemma  2.5,  we  have  that  7t  i— » x(t,  n),  for 
fixed  t  G  I,  is  continuous  and  thus  J  is  continuous  with  respect  to  7C.  We  can,  therefore,  conclude 
that  there  exist  minimizers  for  the  cost  functional  J.  □ 
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Since  the  original  inverse  problem  involves  minimizing  over  infinite  dimensional  sets  Q ), 

pursuing  this  optimization  is  clearly  infeasible  without  some  type  of  finite  dimensional  approxi¬ 
mation.  We  thus  need  to  describe  some  approximation  spaces  over  which  the  optimization  becomes 
computationally  tractable.  Let  Qm  —  {qf}*jL\  be  partitions  of  Q  =  [— r,  0]  for  M  =  1,2, . . .  and 

oo 

Qd  —  U  2m  , 

M=1 

where  the  sequences  are  chosen  such  that  Qd  is  dense  in  Q.  For  a  positive  integer  M ,  let 

{/  M  M  \ 

71  *=  :  n  Pij^qfi  ' '  '  >  52  Pipj^qf  J  > 

M  ) 

qf  e  QM,0<pij,Y,Pij  =  U  =  l,---, > , 

7=1  J 

where  Aq  is  the  Dirac  delta-measure  with  atom  at  q.  That  is,  A q(E)  is  1  if  q  G  E  and  0  if  q  (ji  E. 
Finally,  we  define 

oo 

nD=  Un^. 

M=1 

and  use  the  fact  that  Q  is  a  complete,  separable  metric  space.  By  Theorem  3.1  from  [5],  we 
have  that  Ho  is  dense  in  in  the  pn  metric.  We  can,  therefore,  directly  conclude  that  any 
element  7t  G  Ylc,d  can  be  approximated  by  a  sequence  { 7Lv/; } ,  Kmj  £  such  that  as  Mj  — >  <*>, 
Pn(^My,^)  — >■  0. 

Following  the  discussion  concerning  Theorem  4.1  in  [5],  we  now  state  our  theorem  regarding 
the  continuous  dependence  of  the  inverse  problem  upon  the  given  data  as  well  as  stability  under 
approximation  of  the  parameter  sets  II aci- 

Theorem  2.7.  Let  Q  =  [— r,  0],  assume  that  for  fixed  t  G  [0,7/],  7t  i— >•  x(t,n)  is  continuous  on  Tlaci, 
and  let  Qd  be  a  countable  dense  subset  of  Q  as  defined  above.  Suppose  that  the  obserx’ed  data 
d,„,d  G  M.n,xnd,  nci  <  n,  are  such  that  d,„  — >  d  as  m  —>  °°.  Moreover,  suppose  that  n*M(d„7)  is 
the  set  of  minimizers  for  J(n;dm)  over  K  G  fl^  corresponding  to  the  data  dm.  Similarly,  suppose 
that  n*(d)  is  the  set  of  minimizers  of  J( 7t;d)  over  K  G  Ylad  corresponding  to  the  data  d.  Then, 
ri//(n*M(d,„),rr(d))  ^  0  as  M,m  ->  oo. 
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Proof.  Using  continuous  dependence  of  solutions  on  ti,  compactness  of  Y\ad,  and  density  of  Id/)  in 
Ilad,  the  arguments  follow  precisely  those  for  Theorem  4. 1  in  [5] .  We  thus  can  claim  that 

^(n*M(dw),II*(d))  =  inf  (pn(<,7t) :  <  G  U*M(dm),n  G  LL(d)}  , 
converges  to  zero  as  M.  m  — >  °°.  □ 

Combining  the  results  of  these  two  theorems,  we  can  establish  both  that  there  exists  a  solution  to 
the  inverse  problem  and  that  it  is  continuously  dependent  (in  the  Hausdorff  sense)  upon  the  given 
data.  We  established  well-posedness  of  the  forward  problem  in  §2.2  and  in  this  section,  we  have 
shown  problem  stability  and  a  limited  type  of  method  stability  (stability  under  approximations  of 
the  parameter  sets)  of  our  inverse  problem.  We  can  therefore  conclude  that  our  inverse  problem  is 
well-posed. 

Lastly,  we  note  that  the  theoretical  results  of  this  and  the  next  section  only  apply  to  the  iden¬ 
tification  of  probability  distributions.  The  extension  to  deterministic  parameters,  i.e.,  ones  not 
associated  with  the  delay,  also  follows  readily  from  considering  the  Ordinary  Least  Squares  (OLS) 
optimization  over  J3ad  x  Had  where  is  the  domain  of  admissible  values  for  the  deterministic 
parameters. 


3.  Approximation  Implementation 

3.1.  Approximation  of  the  Forward  Problem.  The  system  described  by  (2.1)-(2.2)  is  a  special 
case  of  a  more  generalized  equation  that  facilitates  a  discussion  regarding  its  approximation.  The 
full  discretization  development,  based  upon  the  ideas  of  [9],  is  presented  in  [6,  18]  and  will  only 
be  summarized  here. 

If  we  denote  a  generalized  solution  space  as  Z  =  M"  x  LniQ'^W'),  we  can  define  the  nonlinear 
operator  srf  :  Qt  (&/)  C  Z  — >  Z  by 

=  {(4>(0),<|>)  G  Z  :  <]>  G  H\Q-,R")} 


^(<K0),4>)  =  (g(<|>(0),<|>),<i>)  • 
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With  this  definition,  we  can  generalize  (2.1)-(2.2)  in  the  form 

z(t)  =  sfz(t) 

(3.1) 

z(0)  =  z  o , 

where  z o  G  Z.  As  we  shall  see  in  §4  below,  the  system  (2.1)-(2.2)  is  a  special  case  of  (3.1)  where 
zo  is  restricted  to  M"  x  PAC(<2,MW). 

With  index  of  approximation  N,  let  {ZN ,  P^,  xZN }  be  our  approximation  scheme  for  (3.1)  satis¬ 
fying  the  conditions  of  Theorem  3.1  in  [9],  where  ZN  is  a  spline  subspace  of  Z,  P'"v  is  the  orthogo¬ 
nal  projection  of  Z  onto  ZN,  and  is  the  approximating  operator  =  P,v.f/P'v.  Thus,  using 
{ Z'v ,  P'v,  .c//N}  we  can  generate  an  approximation  to  the  formulation  described  by  (3.1),  which  we 
denote  by 

ZN(t)  =  £ZNZN(t) 

(3.2) 

zV(o)  =  p NZQ' 

An  alternative  description  of  (3.2)  is 

z?(t)  =Fnzq  +  £  {*/Nz?(s)}ds. 

Theorem  3.1.  Given  the  systems  described  in  (2.1)-(2.2)  and  (3.2)  with  zo  =  (cf>(0).<t>),  e 
PAC(Q;Rn),  under  Condition  2.1,  we  have  that  zN (t)  — >  (x(l),X(),  as  N  — >  uniformly  in  t  on 
the  finite  interval  I,  where  x  is  the  solution  of(2.1)-(2.2). 

Proof.  Following  the  arguments  for  Theorem  2.2  in  [3]  we  readily  obtain  our  desired  conclusion. 

□ 

For  the  full  details,  including  the  selection  of  an  appropriate  basis  for  ZN ,  we  direct  the  interested 
reader  to  [6,  18]  where  piecewise  linear  splines  were  successfully  employed. 

3.2.  Approximation  of  the  Inverse  Problem.  Since  we  do  not  have  direct  access  to  an  exact 
solution  to  (2.1)-(2.2),  our  parameter  estimation  efforts  are  actually  focused  on  solving 

nt 

(3.3)  min  JN(n, d)  =  min  V  \CxN(ti,'K)  —  d,-|“  , 

ftStl aci  TVZTla(i 
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with  nt  data  observations.  As  described  in  the  previous  section,  we  can  obtain  a  uniformly  conver¬ 
gent  (in  t  on  finite  intervals  for  fixed  ft)  numerical  scheme  which  generates  an  approximate  solution 
xN.  We  have  not  yet  shown,  however,  that  as  N  — >  °°,  xN(%j)  — >  .v(ft)  as  7 t/  — >  ft  in  the  pn  metric, 
which  is  needed  for  general  “method  stability”  (see  [5]).  Arguments  similar  to  those  used  in  the 
proof  of  Lemma  2.5  can  be  used  to  show  that  for  fixed  N,  the  approximations  xN  (and  hence  JN  ) 
are  continuous  in  ft  on  IT^.  However,  the  convergence  of  \xN(t1  ft/)  —xN(t,  ft)|  — >  0  as  N,i  — >  °°  is 
not  as  obvious  (but  still  true),  and  is  proven  in  the  following  lemma. 

Lemma  3.2.  Let  tel,  TIE  Llad,  an(J  {ft,-}  e  Ylaci  be  such  that  lim/^oo  pn(ft/,  ft)  =  0.  Then  ifxN(t,Kj) 
is  the  solution  to  (3.2)  under  the  conditions  of  Theorem  3.1,  we  have  \xN(t,Hj)  —  ajV(/.JT)  ^  0  as 
N ,  i  — >  uniformly  in  t  on  1. 

Proof  Denote  P^'  as  the  first  component  of  the  orthogonal  projection  operator  P,v  of  Z  onto  ZN .  By 
using  arguments  similar  to  those  employed  in  deriving  inequality  (2.7),  we  find  that  for  t  €  [0,  tj] 

e(t,x * (ft;),*f  (ft))  =  ^(f, ft,)  —x? (f , ft) |  +  \4 (ft/)  -vf  (ft)L 
<2|vf(ft/)-vf(ft)|oo 

<2  f  Pfg(.VA,(5,ft/),Yf(ft/),ft/)-Pfg(.vyV(5,ft),Yf(ft),ft)Jv 
Jo 

<2  [  |g(//(v,ft/),Yf(ft/),ft/)-g(x/V(5,ft),Yf(ft),ft)|  ds, 

Jo 

and  since  g  satisfies  Condition  2.1,  we  then  know  that 

e(t,j^(ni),j^(n))  <2Kl  [  e(t,4(ni),4  (n))ds  +  ^(ft/,ft) , 

Jo 

where  —  /q  (ft);ft/,ft)ds.  By  the  boundedness  of  |jc^v'(ft) |oo  and  dominated  con¬ 
vergence,  we  have  that  £TN  — >  0  as  i.N  — ►  °o.  A  simple  application  of  Gronwall’s  inequality 

gives  us  that 

e{t,${Tii),$ (ft))  ->0, 

as  iV,  /  — >  and  ft/  — >  ft  in  the  pn  metric,  and  thus  we  conclude  that  \xN(t,llj)  —  jciv(f ,  ft)  |  — >  0  as 


well.  Note  that  the  convergence  is  actually  uniform  in  t  on  I. 


□ 
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Corollary  3.3.  Under  the  conditions  specified  in  Lemma  3.2,  we  conclude  that  \^(t,  tin)  —  x(t,  it)  \  - 
0  as  N  — >  oo  for  TLy  — >  71  in  the  pn  metric.  As  before,  convergence  is  uniform  in  t  on  I. 

Proof.  Consider 

|.XJV(?,7T7v)  -x(t,Tl)\  <  \.^(t,  TIn)  — :TjV(r,7l)|  +  ^  (t  ,Tl)  . 

The  first  term  converges  due  to  Lemma  3.2,  while  the  second  term  converges  as  a  result  of  our 
numerical  scheme  as  given  in  Theorem  3.1.  □ 

With  this  corollary,  we  are  now  prepared  to  examine  questions  concerning  the  existence  of  a 
solution  to  (2.3)  and  (3.3)  as  well  as  the  dependence  of  those  solutions  upon  given  data. 

Theorem  3.4.  There  exists  a  solution  to  the  origincd  and  approximate  inverse  problems  in  (2.3) 
and  (3.3),  respectively.  Moreover,  one  can  find  solutions  to  the  family  of  problems  in  (3.3)  that 
converge  to  a  solution  to  (2.3)  as  N  — >  °°. 

Proof.  As  noted  above,  we  have  that  (naj,pn)  is  compact.  By  Lemmas  2.9  and  3.2,  we  have 
that  both  7t  i— >  x  ( l .  71)  and  71 1— >•  xN  (t,  7t),  for  fixed  tel,  are  continuous  and  thus  both  J  and  JN  are 
continuous  with  respect  to  7t.  We  therefore  know  that  there  exist  minimizers  in  Tlad  for  the  original 
and  approximate  cost  functionals  J  and  JN ,  respectively. 

Let  {7T^}  G  Uad  be  any  sequence  of  solutions  to  (3.3)  and  (tiL  }  a  convergent  (in  pn)  subse¬ 
quence  of  minimizers  (this  is  possible  since  IIaj  is  a  compact  metric  space).  Recall  that  minimiz¬ 
ers  are  not  necessarily  unique,  but  one  can  always  select  a  convergent  subsequence  of  minimizers 
in  n ad-  Denote  the  limit  (in  pn)  of  this  subsequence  as  n*.  By  the  minimizing  properties  of 
k*n  G  T\ud,  we  then  know  that 

(3.4)  JNk(^Nk)  <  JNk(K)  for  all  ji  G  Uad . 

By  Corollary  3.3,  we  have  the  convergence  of  xN(t,llN)  — >  x(t, 7t)  and  thus  JN(n^)  — >  J( 7t)  as 
N  — >  oo  when  pnfTLy.  tt)  — >  0.  Therefore  in  the  limit  as  N, t  — >  °o,  the  inequality  in  (3.4)  becomes 


(3.5) 


7(71*)  <  7(71 )  for  all  7t  G  fl^ , 
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with  7l*  providing  a  (not  necessarily  unique)  minimizer  of  (2.3).  □ 

We  have  proven  not  only  that  there  is  a  solution  to  the  original  and  approximate  inverse  prob¬ 
lems,  but  also  that  as  we  increase  the  state  space  accuracy  of  the  approximate  solution,  in  some 
sense,  it  approaches  a  solution  to  the  original  inverse  problem. 

Following  the  discussion  concerning  Theorem  4.1  in  [5],  we  now  state  our  theorem  regarding 
the  continuous  dependence  of  the  inverse  problem  upon  the  given  data.  We  note  that  the  following 
theorem  describes  general  “method  stability”  in  problems  where  one  approximates  both  the  state 
and  parameter  spaces. 

Theorem  3.5.  Let  Q  =  [— r,  0],  assume  that  for  fixed  t  G  I  =  [0,7/],  7t  i— >•  x(t.7l)  is  continuous  on 
n ad  in  Pn,  XN,  JN  are  the  approximations  given  in  (3.2),  (3.3),  and  let  Qd  be  a  countable  dense 
subset  of  Q  as  defined  above.  Suppose  the  observed  data  d,„,d  G  M."d,  nci  <  n  are  such  thatdm  — >  d 
as  m  — >  °o  and  that  TV$*(dm)  is  the  set  of  minimizers  for  JN( n;dm)  over  7C  G  fl^  corresponding 
to  the  data  dm.  Similarly,  suppose  that  II*  (d)  is  the  set  of  minimizers  of  J( 7C;d)  over  K  G  II aj 
corresponding  to  the  data  d.  Then,  d h ( Tiff  ( dm ) ,  FF  ( d ) )  — ^  0  as  N.M .m  — >  °°. 

Proof.  If  we  combine  the  arguments  of  Theorem  2.7 ,  Theorem  3.4,  and  Corollary  3.3,  as  in  Theo¬ 
rem  4.1  in  [5],  we  readily  can  obtain 

dH(n*NM (dm),n*(d))  =  inf {pn«m,;c)  :  G  Uf(dm),n  G  n*(d)}  -  0, 

as  N,M,m—*°°.  □ 

With  the  results  of  these  two  theorems,  we  can  claim  both  that  there  exists  a  solution  to  the 
inverse  problem  and  that  it  is  continuously  dependent  upon  the  given  data.  We  established  well- 
posedness  of  the  forward  problem  in  §2.2,  and  in  this  section,  we  established  method  stability 
(under  approximation  of  the  state  space  and  parameter  space)  of  our  inverse  problem.  We  can, 
therefore,  conclude  general  well-posedness  of  our  inverse  problem. 
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4.  Example  Illustration 

We  now  outline  an  application  of  this  framework  to  the  HIV  system  (1.1)-(1.4)  of  §1.  For  fixed 
t  El,  r 1  G  M4,  (|)  E  PAC(<2;IR4),  and7t  =  (Pi,  Pi)  £  II  ad,  consider  the  right  side  of  system  ( 1 .  l)-(  1 .4), 

A m  +  Aifi  (4>,  7C)  +  fi (n)  +  f3 (t) , 


where 


and 


— c 

0 

n-c 

0 

0 

nA 

0 

0 

0 

rv  -  3 A 

0 

0 

0 

— Y 

0 

0 

Ai  = 

0 

0 

rv  ~  Sc 

0 

,  M  = 

0 

Y 

0 

0 

0 

0 

0 

ru  st, 

0 

0 

0 

0 

£pi  [<l>2] 

-P1}\1}4 

sPl  [4>2] 

,  fM)  = 

-S  (Lf=t  P\i)  t]2  +pr(iri4 

SPl  [<j)2] 

-S  (Zf=1  r,,)  t!3 

0 

—5  (i:f=i  Tl,-)  Tl4  —  PT1H14 

h(t)  =  [0,0,0,S]r . 

We  note  that  the  nonlinearities  exemplified  by  terms  such  as  pr\iV\A  in  fi  are  both  biologically 
unrealistic  and  fail  to  satisfy  a  global  Lipschitz  condition.  Such  terms  in  ( 1 . l)-(  1 .4)  can  be  correctly 
viewed  as  approximations  to  nonlinear  saturation  terms.  We  address  this  issue,  as  we  did  in  [6], 
with  an  approach  which  is  readily  justified  from  a  biological  viewpoint  by  saturation  dynamics 
principles.  We  define  fi  using  /?  with  standard  saturation  limited  nonlinearities  by  replacing  pv\\V\A 
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by  pi(r|i)r|4  and  8r|;-r|j  by  5,(r|,)r|;-  for  i,j  =  2,3,4,  where  the  functions  p\  and  8;  are  defined  as 


(4.1) 


pi Cni) 


0  ;  rii  <  0 

<  pr |i  ;  0  <  r)i  <  fji 

k  pn i  ;  Tii  <  hi  , 


and 


(4.2) 


Sifti) 


0  ;  Hi  <  0 

<  8h,  ;  0  <  Hi  <  Hi 

k  Shi  ;Hi  <Hm 


(for  finite  upper  bounds  Hi  G  M+,  i=  1 , 2, 3, 4).  The  resulting  function  fa  is  now  globally  Lipschitz 
(see  [3]).  With  g  redefined  for  fixed  t  as 


(4.3)  s(n,<M)  =AH+/i(^,7t)+/2(H)+/3(0) 


the  well-posedness  of  the  inverse  problem  for  it  =  (P\ .  Pi)  follows  if  we  establish  that  g  satisfies 
Condition  2.1.  We  proceed  to  argue  that  Condition  2.1  holds  by  observing  that  for  n,H  G  IR4. 
G  PAC(Q-,R4),  and  n  =  (P1,P2),n=  (Pi, Pi)  e  Uad 

(4.4)  |g(H,^,7t)-g(H,$,7t)|  <  \Ai\ |h-CI  +  |^2|  |/i(^,7t)-/i(^,it)|  +  |/2(H)-/2(H)|  • 


To  bound  the  second  term  on  the  right  side  of  (4.4),  let  us  examine 
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The  last  two  terms  comprise  a  function  which,  for  fixed  (f>2,  converges  to  zero  when  P\ .  P2  con¬ 
verge  to  Pi ,  Pj  in  the  Prohorov  metric,  and  moreover  satisfies  the  dominating  bound  of  Condition 
2.1.  Thus  we  have 


/1  (<l>,  it)  -  /1  (<j>,  *)  |  <  3  |<|>  -  $ L  +  ^(^2;  n,  ft) , 


as  required  in  Condition  2.1. 

To  bound  the  third  term  on  the  right  side  of  (4.4),  we  refer  to  the  arguments  presented  in  Lemma 
4.2.1  in  [6].  The  multidimensional  Mean  Value  Theorem  implies  that  for  T),fj  e  M4 

hiy\)-hiy\)  =  <D/2(ri  +  e(fi-ri)),ri-fi)j0, 

where  the  4x4  matrix  valued  function  is  given  by 


D/2 


^l/2  ^2/2  03/2  ^4/2  5 


where  r),/2  is  the  partial  derivative  of  /2  with  respect  to  the  zth  component  of  its  vector  argument. 
By  the  definition  of  /2  we  know  that  the  quantity  |D/2|  is  be  bounded  by  some  K2 .  A  simple 
application  of  Cauchy-Schwarz  then  yields 


E>/2(ri  +  e  (fi  —  ri))  |  |fi  —  ri I  z/e 


<^2|ri-fi|, 


and  the  combination  of  this  bound  with  (4.4)  yields  our  claim  with  Ki  =  max{3,^2}- 

Therefore,  since  the  function  g  satisfies  Condition  2.1,  we  can  conclude  well-posedness  of  the 
inverse  problem  of  identifying  distributions  associated  with  time  delays  in  the  HIV  viral  infection 
dynamics  described  by  equations  ( 1 . l)-(  1 .4). 
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5.  Conclusions 

Our  efforts  here  are  motivated  by  a  class  of  mathematical  models  which  partition  populations 
into  subclasses  according  to  a  characteristic  temporal  delay.  These  delays  can  be  viewed  as  re¬ 
alizations  of  an  associated  probability  distribution  representing  the  percentage  of  a  subclass  in  a 
population.  The  specific  microscale  structure  of  these  populations  is,  however,  only  accessible 
as  an  expectation  with  respect  to  the  distribution.  In  §1,  we  presented  examples  of  models  ex¬ 
hibiting  these  features  which  occur  in  a  wide  variety  of  fields  such  as  immunological  population 
dynamics,  viscoelasticity  of  polymers  and  rubber,  and  polarization  in  dielectric  materials.  In  §2, 
we  developed  a  mathematical  framework  in  which  we  establish  existence  and  uniqueness  of  the 
forward  problem  and  well-posedness  for  the  inverse  problem  of  estimating  the  probability  mea¬ 
sures.  These  results  included  method  stability  under  numerical  approximations,  thus  leading  to  a 
computationally  feasible  methodology.  Finally,  we  verified  that  one  motivating  model  of  HIV  in¬ 
fection  dynamics  (studied  previously  in  [6])  satisfies  all  the  conditions  of  our  framework,  thereby 
providing  a  theoretical  foundation  for  inverse  problem  computations  with  these  models. 
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